#!/usr/bin/Rscript

# calibrateInflow.R
#
# This scrip preforms the main calibration and validation of the inflow model 
# for the Arcata Waste Water Treatment Plant.  The script is capable of 
# calibrating and valudating two types of models, both of which decompise the
# time series into a seasonal, trend and remainder component. 
# 
#Calibration:
# Model 1 - STL algorithm
#   This model uses the STL algoritm to decompose both the inflow and rainfall 
#   time series and then fit a linear model to the remainder components. For
#   more information on the methosd see:
#       Cleveland, R. B., W. S. Cleveland, J. McRae, and I. Terpenning (1990),
#       STL: A Seasonal-Trend Decomposition Procedure Based on Loess, 
#       Journal of Ofﬁcial Statistic, 6, 3–73.
#
# Model 2 - Fourier model 
#  This model uses fourier analysis to decompose the inflow and time series 

source('callCollect.R')
collect <- 4
resample <- F
usefourier <- F
nsims <- 500
nn <- 25
w <- 14
h <- 7
s <- .9
ciul <- .95
cill <- .05

    #read in and format the data
ifr <- scan('fourier-mean-inflow.txt',quiet=T)
rfr <- scan('fourier-mean-rainfall.txt',quiet=T)
calibFile <- '../overlap/eka-ppt-inch-arc-inflow-mgd-overlap.txt'
predFile  <- '../overlap/eka-ppt-inch-arc-inflow-mgd-overlap.txt'


#---------- data for calibration -----
dat <- as.matrix(read.table(calibFile,fill=T,header=T))
year <- dat[,1]
mon <- dat[,2]
day <- dat[,3]

inflow <- dat[,4]
inflow[ inflow == -9999 ] <- NA
inflow[ inflow == 0 ] <- NA
inflow[ inflow > 10 ] <- NA

rf <- dat[,5]
rf[ rf > 4 ] = NA

    #interpolate NA's (its the best we can do for now)
for(i in 1:nrow(dat)){
    if(is.na(inflow[i])){ 
        inflow[i] <- mean(c(inflow[i-1], inflow[i+1]),na.rm=T)
    }
}
for(i in 1:nrow(dat)){
    if(is.na(rf[i])){
        rf[i] <- mean(c(rf[i-1], rf[i+1]),na.rm=T)
    }
}

    #sum the last collect days
col = callCollect(rf,inflow,collect=collect)
rf = col$x
inflow = col$y

    #create the ts objects
its = ts(inflow, start = c(dat[1,1],dat[1,2]), frequency=365)
rts = ts(rf, start = c(dat[1,1],dat[1,2]), frequency=365)
#------------------

#---------data for validation
dat = as.matrix(read.table(predFile,fill=T,header=T))
year.pred = dat[,1]
mon.pred = dat[,2]
day.pred = dat[,3]

inflow = dat[,4]
inflow[ inflow == -9999 ] = NA
inflow[ inflow == 0 ] = NA
inflow[ inflow > 10 ] = NA

rf = dat[,5]
rf[ rf > 4 ] = NA

    #interpolate NA's (its the best we can do for now)
for(i in 1:nrow(dat)) if(is.na(inflow[i])) inflow[i] = mean(c(inflow[i-1], inflow[i+1]),na.rm=T)
for(i in 1:nrow(dat)) if(is.na(rf[i])) rf[i] = mean(c(rf[i-1], rf[i+1]),na.rm=T)

    #sum the last collect days
col = callCollect(rf,inflow,collect=collect)
rf.pred = col$x
inflow = col$y

    #create the ts objects
its.pred = ts(inflow, start = c(dat[1,1],dat[1,2]), frequency=365)
rts.pred = ts(rf.pred, start = c(dat[1,1],dat[1,2]), frequency=365)
#--------------



year = year[collect:length(year)]
mon <- mon[collect:length(mon)]
wet <- (mon > 11 | mon < 3)
models <- c("STL","Fourier")

rave <- iave <- rper <- iper <- rres <- ires <- list()
a <- b <- fit <- list()
pred <- list(numeric(length(rf.pred)),numeric(length(rf.pred)))
res <- se <- upper <- lower <- list()
z <- qnorm(ciul)

for(j in 1:2){
    cat("\n-----------",models[j],"-----------\n\n")
    if(j==1){  #then stl 
        #decompose ts
        idecomp <- stl(its,'periodic')
        rdecomp <- stl(rts,'periodic')

        rave[[1]] <- mean(rdecomp$time.series[,2])
        iave[[1]] <- mean(idecomp$time.series[,2])
        rper[[1]] <- rdecomp$time.series[year==2000,1]
        iper[[1]] <- idecomp$time.series[year==2000,1]
        rres[[1]] <- rts - rep(rper[[1]],length.out=length(rts)) - rave[[1]]
        ires[[1]] <- its - rep(iper[[1]],length.out=length(its)) - iave[[1]]
    }
    if(j==2){  #then fourier
        rper[[2]] <- rfr
        iper[[2]] <- ifr
        rave[[2]] <- mean(rts-rep(rper[[2]],length.out=length(rts)))
        iave[[2]] <- mean(its-rep(iper[[2]],length.out=length(its)))
        rres[[2]] <- rts - rep(rper[[2]],length.out=length(rts)) - rave[[2]]
        ires[[2]] <- its - rep(iper[[2]],length.out=length(its)) - iave[[2]]
    }



    cat('Remainder Correlation,',models[j],"is:",cor(ires[[j]],rres[[j]]),'\n')
    cat('Remainder Wet Correlation,',models[j],"is:",cor(ires[[j]][wet],rres[[j]][wet]),'\n\n')



    cat('The Average irregular rainfall trend,',models[j],"is:",rave[[j]],'\n')
    cat('The Average irregular inflow trend,',models[j],"is:",iave[[j]],'\n\n')
    
    fit[[j]] <- lm(ires[[j]]~rres[[j]])
    
    cat('The Regression Parameters,',models[j],':\n')
    cat('a =',coef(fit[[j]])[1],'\n')
    cat('b =',coef(fit[[j]])[2],'\n')
    cat("Inflow = a + b * Precip\n\n")
    
    #quartz()
    #plot(idecomp$time.series[,3],rdecomp$time.series[,3],xlab='Inflow Residual',ylab='Precip Residual')
    #abline(fit)
    
    a[[j]] = coef(fit[[j]])[1]
    b[[j]] = coef(fit[[j]])[2]


    theday <- 0
    for(i in 1:length(rf.pred)){
        theday = theday + 1
        if(theday>365) theday <- 1
        pred[[j]][i] = rf.pred[i] - rper[[j]][theday] - rave[[j]]
        pred[[j]][i] = a[[j]] + b[[j]] * pred[[j]][i]
        pred[[j]][i] = pred[[j]][i] + iave[[j]] + iper[[j]][theday]
    }
    res[[j]] <- its.pred - pred[[j]]
    se[[j]] <- sd(res[[j]])
    upper[[j]] <- pred[[j]] + z * se[[j]]
    lower[[j]] <- pred[[j]] - z * se[[j]]

    cat("Standard Error,",models[j],",is:",se[[j]],"\n")

    pred[[j]] <- ts(pred[[j]],start = c(year.pred[1],day.pred[1]), frequency=365)
    upper[[j]] <- ts(upper[[j]],start = c(year.pred[1],day.pred[1]), frequency=365)
    lower[[j]] <- ts(lower[[j]],start = c(year.pred[1],day.pred[1]), frequency=365)


    up = integer(length(pred[[j]]))
    dn = integer(length(pred[[j]]))
    for(i in 1:length(pred[[j]])){
        if(its.pred[i] > upper[[j]][i]) 
                up[i] <- 1
        if(its.pred[i] < lower[[j]][i]) 
                dn[i] <- 1
    }
    cat(round(100*length(up[up==1])/length(up),3),
            '% of observations were above upper CI,',models[j],', expected ',
            100*cill/2,'%\n',sep='')
    cat(round(100*length(dn[dn==1])/length(dn),3),
            '% of observations were below lower CI,',models[j],', expected ',
            100*cill/2,'%\n',sep='')
   cat("--------------------------\n\n") 
}

if(usefourier){
    W <- sapply(se,mean)
    W <- W/sum(W)
    
    Wpred <- W[1]*pred[[1]] + W[2]*pred[[2]]
    Wpred <- ts(Wpred,start = c(year.pred[1],day.pred[1]), frequency=365)
    
    Wres <- its - Wpred
    Wse <- sd(Wres)
    Wupper <- Wpred + z * Wse
    Wlower <- Wpred - z * Wse
    
    cat("Weighted se:",Wse,"\n\n")
}else{
    Wupper <- upper[[1]]
    Wlower <- lower[[1]]
    Wpred <- pred[[1]]
}

postscript(width=7,height=3,file="inflow-ts.eps")
par(mar=c(2.2,2.2,.2,.2))
plot(its,type='l',lwd=2,xlim=c(2005,2008),ylim=c(0,10),xlab="",ylab="")
lines(rts,type='h',col='steelblue')
legend("topright",
        c('Historical Inflow (mgd)','Rainfall (in.)'),
        lty=c('solid','solid'),
        lwd=c(2,1),
        col=c("black","steelblue"))
dev.off()

postscript(width=7,height=5,file="inflow-pred.eps")
par(mar=c(2.2,4.2,.2,.2))
plot(its,type='l',lwd=2,xlim=c(2005,2007),ylim=c(0,10),ylab="Inflow (MGD)")
lines(Wupper,lty='dotted')
lines(Wlower,lty='dotted')
lines(Wpred,col='red',lwd=1)
#lines(rts,type='h',col='steelblue')
legend("topright",
        c('Historical Inflow','Inflow Prediction','95\\% Confidence Interval'),
        lty=c('solid','solid','dotted'),
        lwd=c(2,1,1),
        col=c(1,2,1))
dev.off()


save(rper,rts,rave,iper,its,iave,a,b,z,se,models,file='parms.RData',ascii=T)
cat('Wrote regression parameters and seasonal trend to: parms.RData\n')
save(pred,file='hindcastedInflow.RData')
cat('Wrote hindcasted inflow to: hindcastedInflow.RData\n')
